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Abstract 



A method is proposed to handle the sign problem in the sim- 
ulation of systems having indefinite or complex-valued measures. 
In general, this new approach, which is based on renormalisation 
blocking, is shown to yield statistical errors smaller than the crude 
Monte Carlo method using absolute values of the original mea- 
sures. The improved method is applied to the 2D Ising model 
with temperature generalised to take on complex values. It is 
also adapted to implement Monte Carlo Renormalisation Group 
calculations of the magnetic and thermal critical exponents. 



1 The Monte Carlo sign problem 



In order to evaluate a mult i- dimensional integral 
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using Monte Carlo (MC) one can sample the points in the integration domain 
with a non-uniform distribution, p, which reflects the contribution from the 
measure / at each point, as in the importance sampling This sampling 
gives the following estimate for the integral: 
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The best choice of p is the one that minimises the standard deviation 
squared S. This can be found by variational method leading to the crude 
average-sign MC weight 
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Applying to physical systems, MC method can be used to evaluate the 
expectation value of some measurable quantity 
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Z = J2e~ H , (6) 

where s denotes the dynamical variables and H is the hamiltonian of an 
equilibrium statistical system or the action of some Euclidean quantum field 
theory. 

Since an estimate of the denominator of (^) is independent of a particular 
measurable and is needed in calculating all observables, it is to this that we 
apply the MC method. The integral (H) now assumes the form of the partition 
function Z, upon which the crude weight (^) takes on the explicit form 

Pcrudc - |e-"|/£|e-"|. (7) 

M 

The Boltzmann weight e~ H can in general be real but non-definite, or 
even be complex- valued, in which case it is possible to generalise the absolute 
values of real numbers in the above expressions to those of complex numbers. 
And the error bars now can be visualised as the error radius in the complex 
plane of a circle centred at the complex-valued central MC estimate. The 
variational derivation still goes through as with real numbers. 

We then have in the MC approximation 
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In the above ((sign)) also denotes the average, with respect to a given MC 
weight, of the phase when e~ H is complex. 
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The sign problem H arises when | ((sign) ) | is vanishingly small: then un- 
less a huge number of configurations are MC sampled, the large statistical 
fluctuations of @, because of the small denominator, render the measure- 
ment meaningless. 

Unfortunately, many interesting and important physical problems suffer 
the sign problem like the real-time path integrals of quantum mechanics and 
quantum field theory, lattice QCD at finite temperature and non-zero chem- 
ical potential, lattice chiral gauge theory, quantum statistical system with 
fermions . . . None of the existing proposals is quite satisfactory: complex 
Langevin simulations cannot be shown to converge to the desired distri- 
bution and often fail to do so; others B are either restricted to too small a 
lattice, too complicated, or not general enough or rather speculative. 

In the next section we present another improved method, which is then 
applied to the Ising model in two dimensions and the results will be compared 
with the crude MC of this section, as well as with series-expansion data. 

2 The improved method 

One way of smoothing out the sign problem is to do part of the integral 
analytically, and the remainder using MC ||. The analytical summation is 
not just directly over a subset of the dynamical variables; in general it can 
be a renormalisation group (RG) blocking where coarse-grained variables are 
introduced. We will show below that this does yield certain improvement 
over the crude MC in general. 

Let P{V, V} be the normalised RG weight relating the original variables 
V to the blocked variables V 0, 

P{V',V} > 0, 
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Jp{V',V}dV' = 1. 
Inserting this unity resolution into the integral ([I]) 

dV J dV'P{V',V}f, 

dV'g(V% (9) 

and assuming that the blocking can be done exactly or approximated to a 
good degree such that we then obtain g as a function of blocked variables in 
closed form. An example of the RG blocking which we will employ in the next 
section for the Ising model is the sum over spins on odd sites of the lattice, 
leaving behind a measure g in terms of the other half of the spins on even 
sites. Thus, an MC estimator is only needed for the remaining integration 
over V in (0). As with the crude method of the last section, variational 
minimisation for S of (|3|), with g in place of /, leads to the improved MC 
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This one-step exact RG blocking already improves over the crude average- 
sign method of the last section. Firstly, the improved weight sampling yields 
in (^|) a denominator of magnitude not less than that sampled by the crude 
weight: 
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where we have used the definitions of the sampling weights in the second 
equality, g in (ffi). The inequality is the triangle inequality from the properties 



4 



of P. In other words, from their definitions, ((sign)) improved is proportional to 
((sign)) crude , with the proportionality constant is some function of temper- 
ature and external fields. Both of them vanish when the partition function 
does; away from this point, however, the improved method is no worse than 
the crude one. 

Secondly, it is also not difficult to see that the statistical fluctuations 
associated with improved MC is not more than that of the crude MC, 
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where we have used the definitions of the sampling weights in the second 
equality, definition of g (|9|) in the third. The last inequality is the triangle 
inequality from the properties of P. 

Thus the RG blocking always reduces, sometimes significantly, the statis- 
tical fluctuations of an observable measurement by reducing the magnitude 

of 



| ((sign)) |' 

Note that the special case of equality in ([TTlJl^) occurs iff there was no 
sign problem to begin with. How much improvement one can get out of the 
new MC weight, i.e. how large are the above inequalities, depends on the 
details of the RG blocking and on the original measure /. 
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3 Application to the 2D Ising model 



The Hamiltonian for the Ising model on a square lattice is 

H = -j s n s nl - hJ^Sn- (13) 

(nnf) n 

Here we allow j and h to take on complex values in general. The sum 
over {s} in the partition function is a sum over all possible values of the 
spins s n — {+1, — 1} at site n. The sum over (n'n) is a sum over all nearest 
neighbours on the lattice. For the finite lattice, periodic boundary conditions 
are used. 

The phase boundaries for the complex temperature 2D Ising model with 
h = are found by || 

3 

Re (u) = 1 + 2 a cosu; + 2cos2u; 

3 

Im('u) = 2 a sin a; + 2sin2cu (14) 

where u is taken over the range < u < 2tt, and 

u = e~ 4j . (15) 

In the u plane, this is a limacon, which transforms to the j plane as shown 
in Fig. 1. 

Figure 1: Phase diagram in the complex j plane, h = 0. 
FM=ferromagnetic, PM=paramagnetic, AFM=antiferromagnetic. 



As e~ H j e~ H takes values on the unit circle, the crude MC estimator for the 
denominator of (|8]) might be vanishingly small, but its standard variance S 
is of order unity, leading to the sign problem. 

For the improved method, we adopt a simple RG blocking over the odd 
sites, labeled o. That is, the analytic summation is done over the configu- 
ration space spanned by the o sites; while MC is used to evaluate the sum 
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over the remaining lattice of the • sites. The following diagram shows the 
two sublattices, and how the • sites are to be labelled relative to the o sites, 
for the site labelled x. 



Figure 2: Relative spin positions on a partitioned lattice. 



In general, with finite-range interactions between the spins, one can always 
subdivide the lattice into sublattices, on each of which the spins are indepen- 
dent and thus the partial sum over these spins could be carried out exactly. 
Summing over the spins s Q , 

h s. 



^ — ^ ^ g •sites 

{s,} osites 
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where s+ = sj + + si + s^~. The improved MC weight is then the absolute 
value of the summand on the right hand side of the last expression for Z. 

The quantities to be measured are magnetisation, M, and susceptibility, 
X- These can be expressed in terms of the first and second derivatives of Z 
respectively, evaluated at h = 0. Using the above notation: 



dZ_ 

dh ^ 



•spins 



h ^2 s. r 

e ' sites ]| 2 cosh js^ + h 

osites 



and 



d 2 Z 
dh 2 



]T (sl + tanh(j4 + /*))]] 

.osites J ) 



( r h £ s. 

^2 < e ' sUes f| 2 cosh js+ + h 

•spins (. L osites 

£ (sl + tanh (js+ + h))\ + 

<osites / 



(17) 



x 



E 

osites 



cosh (js+ + h) 



(18) 



7 



4 Numerical results 



In all the simulations, square two-dimensional lattices of various sizes with 
periodic boundary conditions are used. After the RG blocking, half the spins 
go, and the original boundary conditions are maintained. The heat-bath 
algorithm is used to obtained configurations that are distributed with the 
required weights. One heat bath sweep involves visiting every site in the 
lattice once. 



4.1 Autocorrelation 

Two additional benefits arise from the improved method. The first is that 
the number of sites to be visited is halved. While the expressions to be 
calculated at the remaining sites turn out to be far more complicated, the 
use of table look-up means that evaluating them need not be computationally 
more expensive. The second benefit is that correlation between successive 
configurations and hence the number of sweeps required to decorrelate data 
points is reduced. This correlation is quantified in terms of the normalised 
relaxation function (pA (t) for some observed quantity A, 

= ((A(0M( t )»-«A»' 

«1 2 »-«1» 2 1 

The following graph is typical of the behavior near to criticality and 
demonstrates the improvement which is possible. The observable used is the 
real part of the magnetisation versus the number of sweeps. Table 1 shows 
the data used in generating Fig. 3. 

Figure 3: Autocorrelation of the real parts of magnetisation. 
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Table 1: Data for Figure 3. 



Quantity 


Value 


Lattice Size 


32x32 


Total Sweeps 


5,000,000 


Applied Field, h 


+ Oi 


Interaction, j 


0.435 + O.li 


Start 


Cold 


Walk 


Heatbath 



4.2 Improved estimate of << sign >> 

As a test of the improved method, it is compared to the crude one along the 
path OX in Fig. I. 

Figure 4: | << sign >> | vs Re(j), Im(j) = 0.1, h — 

Table 2 shows the data used in generating the remainder of the graphs: 

• Thermalising sweeps is the number of sweeps performed before data is col- 
lected. 

• Data points is the number of configurations used in a measurement. 

• Sweeps between points is the number of sweeps performed between mea- 
surements. 

Note the followings from Fig. 4: 

• Both methods fail close to the phase boundary around Re ( j )= 0.4 
(actually actually no methods can work where Z vanishes) and so results for 
this region are not presented. The important thing is that for a given error, 
the improved method is able to get closer to the phase boundary than the 
crude one. 
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Table 2: Data for Figs. 4-7. 



Quantity 


Value 


Lattice Size 


20x20 


Thermalising Sweeps 


1000 


Data Points 


1000 


Sweeps Between Points 


100 


Applied Field, h 


+ Oi 


Interaction, j 


Re(j) + O.li 


Start 


Cold 


Walk 


Heatbath 



• In agreement with the analytic consideration of section 2, the error bars 
on |<< sign >>| obtained using the improved method are never worse than 
for the crude one. 

• Especially for 0.1 < j < 0.2, the improved method has lifted |<< sign >>| 
drastically, showing that Z ^ but the crude method cannot. This results 
in a big improvement on the satistical errors of the observables, as we will 
see in the following sections. The reason that the value of |<< sign >>| is 
increased more at high temperatures is that the summing of opposite spins 
in cosh [js+ + h] causes a greater reduction in the variance of e~ H / e~ H . 

The gains are more striking if we plot the ratio of the proportional errors, 
R where 

/ 'S'sign 



p _ \«sign»J crudc , x 

Usi ^ n - i s siRn \ ■ 

V «sign» J improved 

This is an important comparison because the errors on the physical observ- 
ables, like magnetisation and susceptibility M and x, depend on the propor- 
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tional error of |<< sign >>|. 

Figure 5: i? sign vs Re(j), Im(j) = 0.1, h = 0. 

4.3 Improved estimate of < M > 

The equivalent graphs for estimates of |< M >| are presented below. For 
clarity, in Fig. 6, only every second data point is shown. The data agrees 
with known behaviour of the magnetisation at high and low temperatures. 

Figure 6: | < M > | vs Re(j) where Im(j) = 0.1, h = 0. 
Figure 7: Ratio of MC error radii of magnetisation. 

4.4 Improved estimate of susceptibility x 

The magnitudes of susceptibility for crude and improved MC are shown in 
Fig. 8. In the region of OX line in the ferromagnetic phase, both methods 
are comparable and consistent with zero. Fig. 9 depicts the ratio of error 
radii of the two simulation methods. 

Figure 8: | < x > I vs R e U) where Im(j) =0.1, h = 0. 
Figure 9: Ratio of MC error radii of susceptibility. 
Comparison with series-expansion data are plotted in Figs. 10 and 11. 
Figure 10: Improved | < x> I vs Re(j). Expansion shown as line. 
Figure 11: Improved < x > vs Re(j). Expansion shown as line. 
The discs in Fig. 11 are the circles of statistical errors for simulated results. 
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4.5 Dependence on lattice size 

Fig. 5 is actually a slice from Fig 12 and Fig 13. These show how the im- 
provement in Rgign depends on the linear lattice dimensions. Two interesting 
trends are apparent: 

• At low temperature the amount of improvement increase with lattice size. 

• At high temperature the amount of improvement decrease with lattice 
size. We do not attempt to explain this behavior, but note that for high 
temperatures one would expect R S i gn to approach some limiting value for 
large lattice sizes. The reason for this is as follows. For both methods, the 
quantity e~ H / e~ H is the sum of arg's over all sites. Hence, by the central 
limit theorem for large lattices it is normally distributed. Its variance is a 
function of the spin statistics, which are not size dependent. 

Figure 12: -R S ign vs Re(j) and linear lattice dimension, 0.1 < Re(j) < 0.3 
Im(j) = o.l, h = 0. 

Figure 13: i? S i gn vs Re(j) and linear lattice dimension, 0.5 < Re(j) < 0.8 
Im(j) = 0.1, h = 0. 



4.6 MC renormalisation group 

We explore the MCRG with both the standard and improved methods at the 
critical temperature on the positive, real axis. It is found that the critical 
exponents of the blocked lattice are the same as those on the original. The 
values of the critical exponents 70 and 71, measured using MCRG are dis- 
played in Table 3. The exact values of 8/15 and 1 are shown at the top of 
the table. 

The data used in generating Table 3 is shown in Table 4: 
• The bootstrap method (JTOj] is used to calculate errors on the critical expo- 
nents. The number of bootstrap samples used, B, is 500. In theory, the limit 
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Table 3: Critical exponents. 





7o (0.533) 




7i (1-00) 




RG iterations 


Crude 


Improved 


Crude 


Improved 


1 


0.532(1) 


0.546(1) 


1.13(4) 


1.08(4) 


2 


0.536(2) 


0.539(3) 


1.12(5) 


1.04(4) 


3 


0.535(4) 


0.539(4) 


1.20(7) 


1.10(5) 


4 


0.535(8) 


0.519(5) 


1.05(7) 


1.11(6) 



of B —y oo should be taken. In practice it is found that the distribution 
changes little for B > 500. 

• The results from the crude and improved methods agree within error. 

• The consistent deviation from the exact value is in agreement with similar 
simulations || and can be explained by truncation of the hamiltonian during 
MCRG and finite size effects. 

• No improvement should be expected (nor is it observed) as there is no sign 
problem in this case. The purpose of these figures is only to demonstrate 
that the improved method is adaptable for use in MCRG. 

Concluding remarks 

We have presented a method towards a partial alleviation of the sign problem; 
it is the earlier proposal in || generalised to include exact RG transforma- 
tions. The sign problem is lessened because of some partial phase cancellation 
among the original indefinite or complex- valued measure after an exact RG 
transformation. 

A particular RG blocking is chosen for our illustrative example of the 
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Table 4: Data for Table 3. 



Quantity 


Value 


Lattice Size 


64x64 


Data Points 


1000 


Sweeps Between Points 


8000 


Applied Field, h 


0.001 + Oi 


Interaction, j 


0.440687 + Oi 


Start 


Cold 


Walk 


Heatbath 


RG Blockings 


5 


Bootstrap Samples 


500 



2D Ising model with complex-valued measure. And this summation over a 
sublattice is the natural choice which always exists for short-ranged interac- 
tions. But other choices of RG blocking are feasible and how effective they 
are depends on the physics of the problems. 

When the quantity to be averaged is not smooth on the length scale of 
the crude weight function, there is an additional source of systematic error 
in the crude, average-sign method. The cancellation in the partial sums may 
reduce this error by reducing the difference in length scales of the measured 
quantities and that of the sampling weights. 
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Table Captions 

Table 1: Data for Figure 2. 

Table 2: Data for Figs. 4-7. 

Table 3: Critical exponents. 

Table 4: Data for Table 3. 
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Figure Captions 



Figure 1: Phase diagram in the complex j plane, h — 0. FM=ferromagnetic, 

PM=paramagnetic, AFM=antiferromagnetic. 

Figure 2: Relative spin positions on a partitioned lattice. 

Figure 3: Autocorrelation of the real parts of magnetisation. 

Figure 4: | << sign >> | vs Re(j), Im(j) = 0.1, h — 

Figure 5: -R S i gn vs Re(j), Im(j) = 0.1, h = 0. 

Figure 6: | < M > | vs Re(j) where Im(j) = 0.1, h — 0. 

Figure 7: Ratio of MC error radii of magnetisation. 

Figure 8: | < % > | vs Re(j) where Im(j) = 0.1, h = 0. 

Figure 9: Ratio of MC error radii of susceptibility. 

Figure 10: Improved | < \ > I vs R e U)- Expansion shown as line. 

Figure 11: Improved < x > vs i?e(j). Expansion shown as line. 

Figure 12: i? sign vs Re(j) and linear lattice dimension, 0.1 < Re(j) < 0.3 

Im(j) = 0.1, h = 0. 

Figure 13: i? S i gn vs Re(j) and linear lattice dimension, 0.5 < Re(j) < 0.8 
Im(j) = 0.1, h = 0. 
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Figure 1 
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